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Abstract 

We apply the Mori-Zwanzig projection operator formalism to obtain an ex- 
pression for the frequency dependent specific heat c{z) of a liquid. By using an 
exact transformation formula due to Lebowitz et al, we derive a relation be- 
tween c{z) and K{t), the autocorrelation function of temperature fluctuations 
in the microcanonical ensemble. This connection thus allows to determine c{z) 
from computer simulations in equilibrium, i.e. without an external pertur- 
bation. By considering the generalization of K{t) to finite wave-vectors, we 
derive an expression to determine the thermal conductivity A from such simu- 
lations. We present the results of extensive computer simulations in which we 
use the derived relations to determine c{z) over eight decades in frequency, 
as well as A. The system investigated is a simple but realistic model for 
amorphous silica. We find that at high frequencies the real part of c{z) has 
the value of an ideal gas. c'{ijj) increases quickly at those frequencies which 
correspond to the vibrational excitations of the system. At low temperatures 
c'{uj) shows a second step. The frequency at which this step is observed is 
comparable to the one at which the a— relaxation peak is observed in the 
intermediate scattering function. Also the temperature dependence of the lo- 
cation of this second step is the same as the one of the a— peak, thus showing 
that these quantities are intimately connected to each other. From c'{lo) we 
estimate the temperature dependence of the vibrational and configurational 
part of the specific heat. We find that the static value of c{z) as well as A are 
in good agreement with experimental data. 
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I. INTRODUCTION 



If a glass-forming liquid is cooled, dynamic observables, such as the viscosity, the diffusion 
constant or the intermediate scattering function, show a dramatic slowing down At 
a certain temperature, the kinetic glass transition temperature Tg, the typical relaxation 
time of the system exceeds the time scale of the experiment and the system falls out of 
equilibrium, i.e. become a glass [Q. In contrast to these dynamical observables, all static 
quantities, such as the volume or the enthalpy, show upon cooling only a relatively weak 
temperature dependence for T > Tg and the glass transition temperature is noticed only by 
a gentle change of slope in their temperature dependence. This change of slope is reflected 
in various susceptibilities, such as the thermal expansion coefficient or the specific heat cy, 
as a relatively sudden drop when the glass transition temperature is reached. The physical 
picture behind this drop in, e.g., cy is that at Tg the structural degrees of freedom fall out 
of equilibrium and the specific heat reduces to the one of a solid having only vibrational 
degrees of freedom. Hence the height of the drop is usually used to estimate that part of the 
specific heat associated with the configurational degrees of freedom. Thus we see that this 
type of experiments can be used to probe the dynamics of the structural degrees of freedom, 
although the information obtained is rather indirect. 

Since the value of Tg is given by the timescale of the experiment, it can be changed if 
the sample is cooled with different cooling rates @J^. Therefore this dependence opens in 
principle the possibility to investigate the temperature dependence of the configurational 
part of the specific heat. However, since the relaxation time of the system changes very 
rapidly with temperature, it is necessary to change the cooling rate significantly (several 
decades) in order to see a significant shift in Tg, which is experimentally rather difficult. 
An alternative way to probe the dynamics of the structural degrees of freedom by means 
of the specific heat is to measure c{z), the (complex) frequency dependence of the specific 
heat of the system in its equilibrium state. Seminal work in this direction was done by 
Nagel and coworkers |6|-p!0| who proposed an experimental setup that allowed to measure 
c{z) also at frequencies as high as 10^ Hz. Independent work in this direction was also done 
by Christensen and others |TI|-|n[. 

Using their measurements of the frequency dependence of the specific heat, Jeong and 
Moon were able to predict cooling and heating curves in differential scanning calorime- 
try (DSC) and found good agreement between their predicted curves and the experimen- 
tal ones |13|. Similar results have been found in molecular dynamics computer simula- 
tions [jl5|,|T6|. Thus these investigations give evidence that the information contained in the 
equilibrium quantity c{z) allows one to understand also the out-of-equilibrium situation that 
is encountered in DSC experiments around Tg. A short review of this can be found in the 



paper by Simon and McKenna [|T7 



The theoretical interpretation of the measurements of c{z) was for some time rather 
controversial [p!^-p3| since the origin of the frequency dependence of c{z) was not really clear. 
Some of these issues could be clarified by experiments done by Birge O and Menon [ITU 



A partial discussion of this dispute is reviewed in Ref. |]T7|. Very recently Nielsen put 
forward a theoretical description of c{z) on purely thermodynamic arguments In the 
present paper we will use a statistical mechanics approach to derive a microscopic expression 
for c{z) and use this expression, which is identical to the one derived independently in 
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Ref. to calculate this quantity from a molecular dynamics computer simulation of 



silica. Interestingly enough this expression has been used before by Grest and Nagel ||25 
to calculate c{z) from a computer simulation of a simple liquid. However, in that paper no 
derivation was given since it seems to have come out "from a stroke of genius" p6|. Due 
to limitations of computer resources the resulting curves for c{z) were unfortunately rather 



noisy, a feature they share with the ones in Ref. [|J]. In contrast to this, great care was 
taken in the present work to obtain reliable data, which in turn allows us to compare the 
temperature dependence of c{z) with the one from other dynamical quantities, such as the 
intermediate scattering function. Thus this permits us to compare the relaxation dynamics 
as probed by c{z) with the one observed by more microscopic methods, such as neutron or 
light scattering experiments. 

The outline of the paper is as follows: In the next section we will use the projection 
operator formalism to derive a connection between c{z) and microscopic variables. In Sec. |T| 
we will give the details of the model and the simulations. In Sec. |V| we will present the 
results and end in Sec. with a summary and discussion. 



II. THEORY 

In this section we derive an expression that relates the frequency dependent specific heat 
to the autocorrelation function of kinetic energy fluctuations. Furthermore we show how the 
thermal conductivity can be calculated from the generalization of this correlation function 
to finite wave-vectors. 

For a dense simple liquid under triple point conditions the only slow modes are the local 
densities of the macroscopically conserved quantities, i.e. the number of particles, the total 
momentum, and the energy. For an A^— particle system these local densities have the form 

N 

X{v,t) = J2Xi{t)S{r-ri{t)), (1) 

1=1 

where for Xi(t) = 1 we obtain the local density fluctuation p(r, t), for Xi(t) = the momen- 
tum density in direction a G {x,y,z}, and for Xi{t) = Ei := /2m+ ^J^jLi ^(ki — rjl) the 
local energy density fluctuation. (Here V{r) is the potential energy.) In the hydrodynamic 
limit the equation of motion for these densities can be derived by using the exact conserva- 
tion laws, some (phenomenological) constitutive equations that relate the spatial derivatives 
of the densities to their currents, and thermodynamic relations between the fluctuations of 
the densities and thermodynamic quantities such as temperature, pressure etc. [0. This 



approach is valid as long as there are no other slow, non-hydrodynamic processes in the sys- 
tem. Well known examples, for which constitutive equations and thermodynamic relations 
have to be modified by, e.g., introducing frequency dependent thermodynamic derivatives, 
are second order phase transitions or systems with internal degrees of freedom ||28|| . 

As mentioned in the Introduction, the dramatic slowing down of the dynamics of liquids 
upon cooling is observed on all length scales. But contrary to second order phase transitions 
no long range order is building up in two particle correlation functions, since the slowing 
down of the dynamics seems not to be caused by a growing correlation length. Instead the 
physical origin for the slow dynamics is the local hindering of the particle motion, the so- 
called cage-effect, i.e. it is due to a mechanism operating on the microscopic scale. However, 
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the consequences of this local slowing down can of course also be detected on mesoscopic 
and macroscopic length scales, e.g. in thermodynamic derivatives, as soon as the time scale 
of the experiment is comparable to the time scale of the structural relaxation. 



In Ref. im a formally exact method to derive frequency dependent constitutive equa- 
tions and thermodynamic derivatives in the canonical ensemble was introduced, which thus 
allows to investigate the consequences of the slowing down of the structural relaxation on 
thermodynamic and hydrodynamic quantities. The main technical tools in that paper are 
thermodynamic response theory and projection operator techniques. In such an approach 
the most important physical ingredient is the appropriate microscopic definition of the tem- 
perature 6{ri, Pi) as a function of the phase space variables rj and pj. From a mathematical 
point of view the choice of a such a function is not unique since the only condition one has 
to fulfill is that 

(0(r„p,)) = r , (2) 

where T is the macroscopic temperature. For the case of supercooled liquids also physical 
considerations impose restrictions to the possible set of functions: In the canonical ensemble 
the temperature of the system is imposed by an external heat bath. Therefore the definition 
of the phase space function 9{ri, pi) should allow that Eq. (H) holds also in the case that T is 
changing on a time scale that is shorter than the one of the structural relaxation. (Such rapid 
changes occur, e.g., in the earlier mentioned experiments on the frequency dependent specific 
heat. But also for the case of a glass below Tg the concept of a temperature is useful, despite 
the fact that the true relaxation time of the system is exceedingly large.) Since structural 
relaxation is related to rearrangements in space, it can be expected that definitions of ^(r^, p^) 
that involve the positions do not fulfill the requirement of being "fast". Momentum, 
instead, can quickly be exchanged between particles, and hence is transfered easily through 
the system. An appropriate microscopic definition of a temperature, obeying the mentioned 
mathematical and physical requirements, can be constructed with the help of the kinetic 
energy of a particle, 

The canonical average {6{pi)) is, as required, the macroscopic temperature T. However, if 
one wants to define a local density related to the temperature, it can of course not be avoided 
that there is a dependence on spatial variables. Thus the best one can do is to make this 
dependence as simple as possible so that it is possible to extract the fast part easily. For a 
classical system a convenient microscopic definition of a local temperature field is therefore 

T(r,t) = $:^(p05(r-r,(t)) . (4) 

i 

The canonical average (T(r, t)) is given by nT, where n = N/V is the number density. Al- 
though it is not a conserved quantity, in supercooled liquids the temperature field T(r, t) has 
a slow component, due to its dependence on the positional degrees of freedom. However, the 
fast fluctuation part can easily be obtained by projecting out these slow density fluctuations: 

5T(r,t)=T(r,t)-Tp(r,t) (5) 
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or in Fourier space 



ST,{t)=T,{t)-Tp,{t) . (6) 

Note that by construction this temperature fluctuation does not have any static couphng to 
functions of the form Gq{r^) = G{r^) exp(iq • Fj), i.e. that depend only on the positions of 
the particles: 

(5r;G,(r^)) = , (7) 

which shows that its dynamic autocorrelation function ^^^(g, t) = (Tg(t)\Tg) cannot be pro- 
portional to any multi-point density correlation function. (For later use we have introduced 
the notation {A\B) := {6A*6B).) Although it is not possible to prove from the Ansatz 
that STq{t) is really a fast variable for g > 0, i.e. that it does not undergo a structural arrest 
in the canonical ensemble, this is something that can be checked in simulations and below 
we will demonstrate that this is indeed the case. Note, however, that the time dependence 
of the temperature fluctuations is of course affected by the slowing down of structural re- 
laxation. This was derived in |2T] within linear response theory. In the following we will 



briefly sketch the main steps in that derivation in order to clarify the physics behind the 
final results. 

Since in this paper we are interested in the frequency dependent specific heat per particle 
at constant volume, the first thing to do is to derive a microscopic expression for this quantity 
and to relate it to the local fluctuations of the temperature. For the sake of simplicity we 
will focus on the case of temperature fluctuations for vanishingly small wave- vectors, since in 
this limit they decouple from the single particle density and current fluctuations. (However, 
in the last part of this section we will also consider finite wave- vectors.) To determine Cv{z) 
we have to calculate, in analogy to the static case AE = c^'^AT, the effect of temperature 
fluctuations on fluctuations of the total energy. For this let us assume that we impose 
adiabatically at time t = a (small) temperature fluctuation ST*{t). In this case the initial 
probability distribution in the canonical ensemble is 



exp{-p{H - h6Tq)) 
Ti exp{-(3{H -h6Tq)) ' 

where h is the conjugate field of 6Tg and /? = l/ksT. 

In the linear response regime the time evolution of any variable Y*{t) is given by 



(8) 



{Y*{t))^E = {Y*it)) + hPiYit)\Tg) , (9) 

where the average on the lefthand side is done in the non-equilibrium ensemble (§), and 
the one on the righthand side are performed in the standard canonical ensemble (/i = in 
Eq. (H)). Using Eq. (|^) for the case Y = Tg we thus obtain for the time dependence of the 
temperature fluctuation 

{6T;{t)) = h P{Tg{t)\Tg) = h P^TTit) , (10) 

and the relaxation of a fluctuation of the energy, Eq. with Xi = Ei, is given by 
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{SE;it)) = h PiE,it)\T,) . (11) 

In analogy to the static case we define the frequency dependent specific heat Cy{z) as the 
coefficient between the fluctuation in the energy and the one in temperature: 

{6E;{z)) = c,{z){6T;{z)) , (12) 

where SEq{z) and STg{z) are the Laplace transforms of SEg(t) and STg(t), respectively. 
Hence we now have to express {Eg{t)\Tq) from Eq. (|TT]) as a functional of 6Tq{t). The 
Laplace transform of {Eg{t)\Tg) is {Eg\R{z)\Tg), with the resolvent R{z) = — z) and 
the Liouville operator C Using a projection operator formalism with the operator P = 
\STq){STq\STg)~^{STq\ to projcct on temperature fluctuations, the correlator {Eq\R{z)\Tq) 
can be written as [1211 



{Eq\R{z)\Tg) = {{Eg\Tg) - {Eg\R' {z)\CTg)} -— {ST^ (z)) . (13) 

XT 

Here R'{z) := Q/{QCQ — z), with Q = 1 — -P, is the reduced resolvent which describes the 
dynamics in the directions of phase space perpendicular to temperature fluctuations (and by 
using the more general formalism from Ref. |2T| also perpendicular to density and current 



fluctuations). The quantity S^rp is defined as 

= (5T;(0)5T,(0)) = 2iVTV3 , (14) 



and is independent of q. (The superscript c indicates the canonical ensemble.) From Eq. (|I3D 
the expression for the specific heat can be read off. We introduce the projection operators 
Pp := \Pq){Pg\Pg){Pg\ sud Q p := 1 — Pp aud use the equations 



CQpE^ = -CQpE^ + 0{q) = -CQE^ (15) 



-QCQ ^ zQ 
z-QCQ ^ z-QCQ ' ^ ^ 

where E^ and Eg are the Fourier transform of the kinetic and potential energy fields E^{v) 
and E^{t), respectively [|29|. The specific heat can now be expressed as 

c.(^)=C + ^lim^«|i?'(^)lO , (17) 

where c^'^ = Ct, = ■^\img^Q{Eg\Qp\Eg) is the equilibrium value of the specific heat per 
particle. We emphasize, that Eq. ([17D is an exact expression for the frequency dependent 
specific heat at constant volume. It only relies on Eq. (^ as a physical reasonable definition 
of temperature fluctuations. 

Expression (|T^ reveals the physical origin of the frequency dependence of the specific 
heat very clearly: During the process of structural relaxation the system probes the potential 
energy landscape. This dynamics is governed by the (slow) relaxation of the densities and 
their higher order products and hence gives rise to a nontrivial slow dynamics even in the 
projected dynamics R'{z), although the latter does not contain any hydrodynamic poles. 
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Expression ([TTD is the ideal starting point for approximation schemes hke the mode 
couphng theory. However, for exact calculations of Cy{z), e.g. in computer simulations, it 
cannot be used, since it is not possible to implement the projected dynamics. To solve this 
problem we will now express Cy{z) in terms of correlation functions that can be measured 
in a real dynamics, such as ^TT{(l,t). For this we make use of Refs. pT] , pT| where similar 
calculation as just presented for Cv{z) were done for Cp{z), the specific heat at constant 
pressure, as well as the frequency dependent heat conduction coefficient X{z). Using these 
results it can be shown that for small wave- vectors the Laplace transform of t) obeys 

the exact equation of motion 



^TT{q,z) :=z dtexp{tzt){5T*{t)5T,iO)) 



zcp[z) + \[z)q - z^ 



(18) 



z"^ — q^KB{z)/mn 



where Kb{z) is the frequency dependent bulk modulus and n is the particle density pT| , |3T . 



Although in general A depends on frequency z, theoretical arguments show that in su- 
percooled liquids this dependence is only weak, in agreement with experimental findings ^ . 
Therefore X{z) can be replaced by its static value, X{z = 0) = iX/n, where A is the equilib- 
rium heat conduction coefficient ||^. Note that the form for the equation for ^tt{(1,z) is 
exactly as the one in linearized hydrodynamics |^ except that now the transport coefficients 



and the thermodynamic derivatives are generalized to functions of the complex frequency z. 

In the limit of vanishing wave- vectors we thus obtain from Eq. (|T^) the following relation 
between ^tt{<1,z) and Cv{z): 

lim5^#^ = ^^ . (19) 



Sq^rp 2ZC„(Z) 



We emphasize that it occurs in Eqs. (|l^) and (|l^) is exactly the dynamic specific 

heat from Eq. (0) that relates fluctuations in temperatures to the ones in energy. Using 
Eq. (0), we thus can express this dynamic specific heat by a standard correlation function: 

cy{z)= -^^f'^TT (20) 

Therefore we find that in the canonical ensemble the specific heat can be expressed by 
the autocorrelation function of the temperature fiuctuations, a quantity which is readily 
accessible in simulations. (An equivalent definition of the specific heat in terms of potential 
energy fiuctuations is given in ||33[|). 



It is important to note that the expression (|T8D for (^tt is valid only in the canonical 
ensemble. To obtain the corresponding relation in the microcanonical ensemble, a general 
transformation due to Lebowitz et al. can be used [^, which relates averages of time depen- 
dent functions in the canonical ensemble to their averages in the microcanonical ensemble. 
Denoting 5T{q = 0) by 6To we thus obtain: 



Mcv 



21 r^.rr.. 2 



d{To) 



dp 
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Here we used, that (Tq) = A^T [35]. If we set t = 0, the equihbrium specific heat c^'^ can 
thus be written as 



3k 



B 



e, 1-^(0) , (22) 



2d 

where we have introduced the normahzed temperature autocorrelation function 

K{t):=^^TT\t)ISTT ■ (23) 
Putting this back into Eq. ( pT]) and using Eq. (|20|) we finally obtain 

^-(^) = 1 _ ;r(t = 0) -^LT[ir(t)](z) ' ^^^^ 

where LT stands for the Laplace transform. 

Equation (0), with a slightly different definition of the function K{t)^ was first used in 
||25| as an ad hoc generalization of Eq. (^) to non-vanishing frequencies and later derived 
within thermodynamic response theory by Nielsen 



So far we have focused on the autocorrelation function of Tq in the limit g — > 0. In the 
following we will now consider the case g > 0. From Eq. ([T8|) it follows immediately that 
for small but not vanishing q the decay of <^tt{<1i t) is dominated by the hydrodynamic pole 
at uj = —iDTq^ with Df = \/{nc'^). Below we will see that it is important to consider also 
corrections of order to Dt- As it will turn out, the frequency dependence of the specific 
heat will give the main contribution in these corrections. From theoretical considerations 
pl| as well as from experiments ||^ we know that the heat conduction coefficients do not (or 
only very weakly) depend on frequency. In analogy to the specific heat at constant volume, 
Cv{z), the specific heat at constant pressure, Cp(z), will also depend on frequency. Using 
a generalization of Legendre transformations between different thermodynamic ensembles 
it can be shown that Cp{z) can be written in the form Cp{z) = d:^ + zA(z), where A(z) 
is a frequency dependent function with a positive spectrum ||31[]. For 2; = u; + ze, it is 
given by Xmi^^^l^ioj + ie) = A'{u) + i/S"{uj). The value of A"(a;) at = has several 
contributions, one being proportional to the longitudinal viscosity rji = rjv + ^rjs, where r]v 
and rjs are the bulk and shear viscosity, respectively. The relaxation time Tq of $tt(5', ^) 
which is related to the hydrodynamic pole can be calculated by determining the frequency 
at which the main denominator of Eq. (|T8|) vanishes up to order q^. Making the Ansatz 
1/Tq = ujr = —iDxq^ + iBq'^ — Aq^ we obtain 



^ 2 / ^"{^ = 0) ^ 2 7 - 1 ..2 2\ ..2 ^'{^ = 0) 4 



where 7 is Cp^/c^'', 7 — 1 is the Landau - Placzek ratio, and Kb is the static bulk modulus. 
The real part indicates an experimentally undetectable shift of the position of the Rayleigh 
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peak to finite frequencies. Since it is tlie frequency of an overdamped harmonic oscillation, 
it does not influence the relaxation time Tq of the correlator in real time. From Eq. (^Sj) the 
relaxation time up to order of the exponentially decaying temperature autocorrelation 
function is thus given by 

1 1 A"(cj = 0) ^ 7 - 1 2x 

ojR A Cp KBj/mn 

It is interesting to note that for the case of supercooled liquids the term proportional 
to q^ is only positive due to the existence of a frequency dependent specific heat, i.e. that 
A" > 0. Without it, the subtraction of the positive Landau-Placzek ratio 7 — 1 would 
lead to a negative offset in a plot of Tg against Below we will check the validity 

of this g— dependence and calculate from this relation the value of A. Note that in the 
case of finite wave-vectors the autocorrelation functions {6T*{t)6Tq{0)) in the canonical and 
microcanonical ensemble are identical, and hence it is not necessary to use the transformation 



formula of Lebowitz et al pl|. 



III. MODEL AND DETAILS OF THE SIMULATIONS 

In this section we discuss the system we used to test the ideas presented in the previous 
section and give the details of the simulations. At the end we also briefly discuss the influence 
of finite size effects on the results. 

Although the formalism presented in the previous section is of course applicable to all 
types of glass forming liquids, it is of special interest to investigate a system which exists 
also in reality, since this opens the possibility to compare the results from the simulations 
with those from experiments. To do a simulation of a real material one needs a potential 
that describes reliably the interactions between the atoms of this substance. In the case of 
silica such a potential does indeed exist, since a decade ago van Beest et al. (BKS) used ab 
initio calculations to obtain a classical force field for this material |26|. The functional form 
of this potential is given by 

0^^(^) = M^ + A,^exp(-5,^r)-^ «,/?e[Si,0], (27) 

where r is the distance between two ions of type a and [3. The value of the constants 
qa,q/3, Aa/3, Ba/3, and Ca/s can be found in Ref. The short range part of the potential 
was truncated and shifted at a distance of 5.5A, which leads to a better agreement of the 
results for the density of vitreous silica as predicted from this model with the experimental 
values 0. In the past it has been shown that this potential is able to reproduce reliably 
various properties of amorphous silica, such as its structure, its vibrational and relaxational 
dynamics, the static specific heat below the glass transition temperature, and the conduction 
of heat |]5|j37|-^ . Thus it is reasonable to assume that this potential will also reproduce well 



the quantities needed to calculate the frequency dependent specific heat, i.e. the correlation 
function K{t) in Eq. (|23|) . 

Using the BKS potential, the equations of motion were integrated with the velocity form 
of the Verlet algorithm with a time step of 1.6fs. The sample was first equilibrated by cou- 
pling it every 50 time steps to a stochastic heat bath for a time which allowed the system 
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to relax at the temperature of interest. (We assumed a sample to be relaxed if the co- 
herent intermediate scattering function at the wave-vector corresponding to the main peak 
in the static structure factor had decayed to zero within the time span of the equilibra- 
tion 1^.) After this equilibration we started a microcanonical run for the production from 
which we determined the equilibrium dynamics at various temperatures. The temperatures 
investigated were 6100K, 4700K, 4000K, 3580K, 3250K, 3000K, and 2750K. At the lowest 
temperature K{t) decays so slowly that runs with 30 million time steps were needed to 
equilibrate the sample, which corresponds to a real time of 49ns. All the simulations were 
done at a constant density of 2.36g/cm'^ which corresponds to a pressure around 0.87 GPa 
T5| . Since the temperature expansion coefficient of silica is small PJ^, it can be expected 



that a constant pressure simulation would basically give the same results as our simulation. 

To obtain also results in the glass at room temperature we cooled the sample from 
3000K to 2750K with a cooling rate around 3 ■ lO^^K/s, and then rapidly quenched it to 
300K. The so obtained glass was annealed for additional 300,000 time steps before we started 
the measurement of the various quantities. 

Since the main observable of interest, K{t), is a collective variable, it has a relatively 
large statistical error. Therefore we averaged all the runs over 200 independent samples for 
T > 4700K, over 100 samples for 4000K > T > 3000K, and 50 samples for T = 2750K. In 
order to make so many independent runs the system size had to be rather small. Therefore 
we choose a system of 112 silicon and 224 oxygen ions in a cubic box of size 18. 8A. The 
Coulombic part of the potential was evaluated with the Ewald method and a parameter 
a of 7.5A~^. All these calculations took around 8 years of single CPU time on a parallel 
computer with high end workstation processors. 



As it was shown in Refs. |^ , |50| , the dynamics of network glass formers shows quite strong 
finite size effects since the small systems lack the acoustic modes at small wave-length. Thus 
it can be expected that if we determine K{t) for a system of 336 particles the result will be 
different from the ones for a system of macroscopic size. In order to check the influence of 
the system size on our results we have made some test runs with 1008 particles and found. 



in agreement with the results of Refs. [47,50|, that at long times the dynamics of the larger 



system is a bit faster than the one of the small one However, the qualitative behavior 
of the various relaxation functions are independent of the system size and therefore we can 
expect that the results presented in this work will hold also for larger systems. 



IV. RESULTS 

In this section we present the results from our simulation, i.e. the temperature depen- 
dence of K(t) and the frequency dependence of the speciflc heat. At the end we briefly 
discuss the time and temperature dependence of the generalization of K(t) to flnite wave- 
vectors. 

As it is obvious from Eq. (|2^) , the flrst step in the calculation of the frequency dependent 
speciflc heat is to determine K{t), the autocorrelation function of the fluctuation in the 
kinetic temperature. Since K{t) is a collective quantity it is relatively difficult of determine 
it with high accuracy. This is shown in Fig. [I|a where we plot K{t)/K{0) at T = 3000K. 
Each of the thin curves corresponds to the average over ten independent samples. From the 
figure we recognize that even if the shape of the different curves is quite similar, their height 
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varies considerably. Thus one should realize that even if the average over 100 samples gives 
a quite nice and smooth curve (bold solid curve) it still might be subject to a significant 
statistical error. 

In Fig. [p3 we show the time dependence of K(t)/K{0) for all temperatures investigated. 
From this figure we see that for all temperatures this correlator decays within lO^^ps to a 
value smaller than 0.2, i.e. very rapidly. For the highest temperatures K{t) then shows a 
small shoulder at around O.Olps, a feature which, at these temperatures, is not observed 
in correlation functions such as the intermediate scattering function F{q,t), where q is the 
wave-vector p2| , ^ . (See Ref. for a definition of F{q,t)). With decreasing temperature 
this shoulder extends to larger and larger times until we observe at the lowest temperature 
a well defined plateau which extends over several decades in time. Thus from this point of 
view K(t) behaves qualitatively similar to F{q,t). Since for this quantity it is customary 
to refer to this final decay as the "a— relaxation" we will use the same term in the case of 
K{t) as well. Note that the height of the plateau in K{t) is only around 0.1, which shows 
that most of the correlation of the kinetic energy is lost during the brief ballistic fiight of 
the particles, t < 10~^ps, inside the cage. 

From the figure we also see that with decreasing temperature the correlation function 
starts to show a minimum at short times. The reason for this feature is that at low temper- 
atures and short times the system behaves similar to a harmonic solid and thus it can be 
expected that the connection between K{t) and the velocity-autocorrelation function J{t), 
see Eq. ( [A7|) in the Appendix, starts to hold. (We remind the reader that in supercooled 
liquids the velocity-autocorrelation function shows a dip at short times.) That this is in- 
deed the case is demonstrated in Fig. |I]c, where we compare the two correlators at a high 
and a low temperature. We see that for T = 300K, i.e. in the glass where the harmonic 
approximation is valid, K{t) and J(2t) are identical within the accuracy of the data (solid 
lines). For temperatures at which the system is still able to relax the situation is different. 
At short times K{t) (bold dashed line) shows oscillations with extrema which are located at 
the same times at which also J(2t) (thin dashed line) shows maxima and minima. Thus the 
harmonic-like character of the motion on these time scales is clearly seen. For larger times, 
however, J{2t) goes rapidly to zero whereas K{t) shows the above discussed plateau before 
it decays to zero at very long times. 

In order to investigate this point in more detail we can make use of Eq. ( [A8| ) , which 
relates the spectrum of K{t), K{uj), to g{oj), the time Fourier transform of the velocity- 
autocorrelation function. At low temperatures the latter quantity is nothing else than the 
density of states of the system. (In order to calculate these Fourier transforms we made 
use of the Wiener- Khinchin theorem which relates the power spectrum of a time dependent 
signal to the Fourier transform of the corresponding autocorrelation function ||51]].) In Fig. |^ 
we show K{uj) for a temperature in the melt and in the glass (dashed lines) and compare 
these curves with g{uj/2)/8 at the same temperatures (solid fines). We see that, within 
the accuracy of our data, in the glass the two curves are indeed identical. (We remind the 
reader that the two peaks at high frequencies correspond to intra-tetrahedral motion of the 
atoms, whereas the broad band at lower frequencies corresponds to (mostly) delocalized 



inter-tetrahedral motion i53|J531.) For the case of the melt, however, the two curves differ 



significantly from each other. Although the general shape of K{uj) and g{oo/2) are similar, 
the former quantity has a much smaller intensity at the two peaks at high frequencies but a 
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higher one in the broad band. Note that this "excess" in low frequency modes has nothing 
to do with the fact that K{t) shows at long times a a— relaxation, whereas J{2t) does not, 
since at T = 3000K the typical frequencies of the a— relaxation are smaller than ITHz, 
and their contribution to the spectra can only be seen as the narrow peak at = (see 
Fig. 0) [^. Thus from this figure we can conclude that at low temperatures the harmonic 
approximation is very good whereas at intermediate temperatures significant deviations are 
observed. 

Since c{z) is related to the Laplace transform of K{t), see Eq. (p4D, one has to calculate 
this transform with high accuracy. From Fig. ^3 we see however, that at low temperatures 
K{t) extends over many decades in time which makes the calculation of this transform a 
non-trivial matter. From this figure one also recognizes that, despite the large number of 
samples we used, the curves have still a significant amount of noise, most noticeable at long 
times and the lowest temperatures, since for these we used fewer samples. Since within 
the accuracy of the data the shape of the curves does, in the late a— relaxation regime, not 
depend on temperature, we substituted for T < 4700K that part of the curve which was 
below 0.02 by the corresponding part of the curve for T = 6100K, after having it shifted to 
such large times that the resulting curve was smooth at 0.02. Subsequently the so obtained 
curves were smoothed and the Laplace transform calculated making use of the formula by 
Filon [g. 

In Fig. |]we show the real and imaginary part of c{z) for all temperatures investigated. (In 
the following we will assume z = u + ie, with e — 0.) To discuss the frequency dependence of 
these curves let us focus for the moment on c'{uj) for T = 2750K, the lowest temperature at 
which we could equilibrate our system. At very high frequencies c'{uj) becomes independent 
of LJ and has a value of 1.5, which is the specific heat of an ideal gas. This result is reasonable 
since these high frequencies correspond to times at which the dynamics of each particle is not 
affected by the other ones, i.e., they move just ballistically. Upon lowering the frequency the 
specific heat rises rapidly since we are now in the frequency regime in which the dynamics of 
the system is mainly dominated by vibrations, see Fig. Therefore at these frequencies the 
system can take up energy giving rise to the increase of c'{uj). Note that before this increase 
occurs, c'{uj) shows a little dip around 40THz, i.e., it falls below the ideal gas value. This 
dip can easily be understood by the harmonic picture proposed above, because the specific 
heat for a single harmonic oscillator shows a singularity at its resonance frequency. Since in 
our system we have many different oscillators that have typical frequencies between 1 and 
80THz, this singularity is smeared out and results in the dip and subsequent strong increase 
of c'{uj). The same mechanism is the reason for the little peak at around 5THz. 

If the frequency is decreased further, c'{uj) stays constant until uj~^ is on the order of 
the time scale of the a— relaxation, i.e. the time scale of the structural relaxation. In 
this frequency range the system is again able to take up energy and hence c'{oo) increases 
again. At even lower frequencies c'{uj) becomes constant, i.e., it has reached the value of 
the static specific heat. This sequence of features in c'{u) can of course also be found in 
c"{uj), since the two functions are related by the Kramers- Kronig relation. In Fig. |^ we see 
that at high frequencies the imaginary part has a peak which corresponds to the vibrational 
excitations in the system. At much lower frequencies we find the a— peak which corresponds 
to the structural relaxation process, i.e. the type of dynamics in which particles change their 
neighbors. For future reference we will introduce the terms "vibrational and configurational 
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part of the specific heat" by which we mean the height of the plateau in c'{uj) at intermediate 
frequencies and the height of the step at low frequencies, respectively, and will denote them 
by Cconf and cvib- 

Let us now discuss the temperature dependence of c'iui) and d'iuj). From Fig. ^ we see 
that the specific heat at intermediate frequencies is essentially independent of T, since the 
vibrational motion of the ions is just a weak function of temperature. The main effect of an 
increase in temperature is that the height of the flat region at very small cu, i.e. the static 
specific heat, increases, and that the crossover from this region to the plateau at intermediate 
frequencies moves to higher frequencies. At the highest temperature this crossover frequency 
has moved up to such high frequencies that no intermediate plateau is visible anymore, which 
means that in the system there is no longer a separation of time scales for the vibrational 
and relaxational processes. Furthermore we find that at the highest temperature the height 
of the plateau at small frequencies has decreased, i.e. that the value of the static specific 
heat has decreased. This effect is most likely related to the fact that silica shows a density 
anomaly, which for the present model occurs at around 4600K |^. 

The temperature dependence just discussed is also found in the imaginary part of c{z) in 
that, with increasing temperature, the a— peak moves continously to higher frequencies until 
it merges with the microscopic peak. All this behavior is qualitatively similar to the one 
found for dynamic observables that measure the structural relaxation, such as the dynamic 
susceptibility |]57| , |58| . Thus this gives evidence that the observable related to the structural 
relaxation and c{z) are closely connected to each other. 

We also note that at low temperatures the form of the curves at low frequencies as well as 
their temperature dependence resembles very much the ones found in experiments |6|- p!0| , p!3 



The main difference is that in the simulation it is possible to measure c{z) even at such high 
frequencies that the effect of the microscopic vibrations becomes visible. Thus it is possible 
to follow continuously the evolution of c{z) from the microscopic regime to the mesoscopic 
one, i.e. to investigate the whole frequency dependence of c{z) from the liquid state to the 
viscous state. In contrast to this, experiments can probe only the frequency range below 
lO^Hz and therefore only the a— regime is observable. However, since experimentally it is 
much simpler to equilibrate the material also at temperatures close to the glass transition 
temperature, c{z) can be measured at significantly lower temperatures than in a simulation. 

Finally we mention that we have included in Fig. |^ also the data for T = 300K, i.e. a 
temperature at which the system is deep in the glass state. We see that this curve follows the 
pattern of the equilibrium curves very well in that it shows also the "harmonic resonance" at 
high frequencies and then a plateau at lower frequencies. No second plateau is seen in c'{uj) 
(or an a— peak in c"{uj)) at very small frequencies since at this temperature these features 
would occur at such low uj that they are not visible within the time span of the simulation 
(or even an experiment). 

Since the part of the specific heat that is related to the structural relaxation is the 
a— relaxation peak at low frequencies, we will study this peak now a bit in more detail. In 
Fig. ^ we show an enlargement of this peak for intermediate and low temperatures. We 
clearly see that with decreasing temperature the area under the peak becomes smaller, 
which means that the configurational part of the specific heat decreases. This temperature 
dependence might be somewhat unexpected since for other dynamic susceptibilities, such as 
the one connected to the intermediate scattering function, one finds that the so-called time 
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temperature superposition principle is valid, i.e. a decrease in temperature just gives rise 



to a horizontal shift of the a— peak |^ , |58[| . However, since the area under the a— peak is 
related via the Kramers-Kronig relation to the height of the step at low frequency in c'(a;), 
i.e. the configurational part of the specific heat, it is not surprising that this area depends on 
temperature, and below we will show that this is indeed the case. This tendency can also be 
easily understood from Eq. (^l]) : For this we assume that in the a— relaxation regime at low 
temperatures the function K{t) can be written as K{t,T) = f(T)K(t/T(T)), where /(T) 
is the height of the plateau and t(T) is the typical relaxation time. That this assumption 
is reasonable can be inferred from Fig. |I]b. It is now simple to show that c"{u!) can be 
approximated by c"{uj) ^ c"(cur(T))/(T)/(l - K{t = 0)f = g"(cjr(T))/(T)4(c^^)V9fc|, 
with a master function c". (Here we also made use of Eq. (|2^) ). Since the function c" is 
assumed to be independent of temperature, we see that the whole T— dependence of c"{u) is 
given by a shift in frequency proportional to t(T) and a vertical rescaling by f{T){cl'^y . Thus 
we conclude that the configurational part of the specific heat is proportional to f{T){d^^Y. 
For the present system both, f{T) and c^'^, decrease with decreasing temperature, and 
thus it is clear that the same is true for Cconf- However, in certain materials, such as 
fragile glass-formers, it is sometimes observed that the specific heat increases with decreasing 
temperature. Thus in these cases the temperature dependence of Cconf does not have to 
decrease monotonically, but it might, e.g., exhibit a local maximum. 

Although the height of the a— peak changes, its shape seems to be independent of tem- 
perature. To demonstrate this we plot in the inset of Fig. ^ c'' {u) / c'' {urnax) versus cj/cUmax, 
where uj^a_y^{T) is the location of the a— peak. Since the curves for the different temperatures 
fall on top of each other, to within the accuracy of the data, we conclude that the shape does 
indeed not change with temperature. Also included in the figure is the Fourier transform of a 
Kohlrausch- Williams- Watts-law with a stretching exponent 0.9. We see that this functional 
form fits the master curve quite well, at least if one does not go to too high frequencies. At 
these higher frequencies the scaling breaks down due to the presence of the microscopic peak. 
We also mention that the (relatively large) value of the stretching parameter is reasonable, 
since in strong glass formers the stretching in the structural relaxation is usually weak, and 
indeed we have found that also for the present model the structural relaxation shows only a 



weak stretching 42 



Further evidence that the structural relaxation and the frequency dependent specific heat 
are closely connected to each other can be obtained by comparing the typical time scales 
for these functions. For this we have calculated the (incoherent) intermediate scattering 
function Fs{q,t) |^ for a wave-vector q = 1.7A~^, which corresponds to the location of 
the first peak in the static structure factor [^. We have defined the a— relaxation time 
Tp, a G {Si, O}, by the time it takes this correlation function to decay to 1/e of its initial 
value. To characterize the time scale for the specific heat we have determined from c"{uj) 
the position of the maximum of the a— peak and defined the relaxation time Tc as l/umax- 
In Fig. 1^ we show the temperature dependence of these relaxation times in an Arrhenius 
plot. From that graph we see that the relaxation times Tc and Tp track each other very 
well in that at low temperatures both of them show an Arrhenius law with a very similar 
activation energy (numbers are given in the figure). With increasing temperature, deviations 
from this law are seen, the origin of which have been discussed elsewhere , but also these 
deviations are the same for both quantities. That the structural relaxation and the specific 
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heat do indeed track each other is demonstrated in the inset, where we plot the ratios 
TpjTf. as a function of inverse temperature. Since we see that these ratios are independent 
of temperature to within the accuracy of the data, we can conclude that the temperature 
dependence of the three quantities is indeed the same. 

In the discussion of Fig. |^ we have mentioned that the increase in c![lS) at high fre- 
quencies is due to the vibrational degrees of freedom, whereas the step at lower frequencies 
is due to the relaxation of the configurational degrees of freedom. By measuring the heights 
of these two steps it thus becomes possible to determine the contribution of the vibrational 
and configurational degrees of freedom to the static specific heat and to investigate how 
these quantities depend on temperature |^. The results obtained are shown in Fig. |^, 
where we plot c^'', the value of diuS) at very low frequencies which is hence the static spe- 
cific heat, Cconf, the height of the step in d{uS) at low frequencies, and Cvib, the value of 
d{uS) at intermediate frequencies. Several observations can be made: Firstly, Cvib shows a 
very regular temperature dependence which can be approximated well by a linear function 
of temperature, at least in the temperature range investigated. An extrapolation of this 
temperature dependence to lower temperatures shows that Cvib attains the harmonic value 
of only at low temperatures (~ lOOOK), i.e. at temperatures that are well below the 
(experimental) glass transition temperature, which is at 1450K a value that seems to 



be reproduced reasonably well by the present model [^]. Hence we find that Cvib is affected 
by anharmonic effects even at relatively low temperatures. We also mention that the vibra- 
tional part Cvib in the temperature range 2750K <T< 3500K seems to be in nice agreement 
with a linear extrapolation of the experimental specific heat of the glass below Tg = 1450K 
to higher temperatures, i.e. if one leaves out the increase of the specific heat due to the glass 
transition. 

The temperature dependence of Cconf is much more pronounced in than the one of Cvib, in 
that it shows around 4000K a crossover from a relatively weak temperature dependence at 
high T to a stronger one at low T. Furthermore we see that Cconf is significantly smaller than 
Cvib, which is in agreement with the experimental result that strong glass formers show only 
a small drop in the specific heat when the temperature is lowered below the glass transition 
temperature, i.e. when the relaxational degrees of freedom are frozen ||63|. We also note that 
for a strong glass former one expects that the Kauzmann temperature is very small []64 |. 
From the figure it seems however, that a naive extrapolation of Cconf to lower temperatures 
leads to a intercept with the temperature axis around Tconf ~ 2000K. Since the inequality 
T'conf < Tk must hold, this type of extrapolation thus leads to an estimate of > 2000K. 
This high estimate of is corroborated by recent results of the same model in which Tk 
was estimated by the direct calculation of the entropy and a subsequent extrapolation to 
lower temperatures |l65| . 

These results for Tk depend of course crucially on the way Cconf is extrapolated to lower 
temperatures. From Fig. ^ it is clear that it is also possible to make this extrapolation in such 
a way that, e.g., at T = 1450K its value is around O.S/cs/particle, i.e. equal to the height 
of the step in the experimental curve at Tg (see experimental curves in figure) |^|. If the 



extrapolation is done in this way, the estimate of Tconf is moved to much lower temperatures. 
Thus it will be very interesting to attempt to do simulations at even lower temperatures in 
order to minimize the effects of this extrapolation. For this it will of course be necessary to 
equilibrate the system at even lower temperatures, which is computationally difficult. One 
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promising way to achieve this is the so-called method of "parallel tempering" , and work in 
this direction is presently done [l67|-|69| . 



Also included in the figure is the specific heat of the system as calculated from the 
harmonic approximation |^ (dashed line). This was done by determining the eigenvalues 
of the dynamic matrix, and hence the density of states g{uj), and using the expression 

_^ .00 .^,(.)exp(/.c./fc.T)^^_ 



More details on this calculation can be found in Ref. where it has been shown that 
this theoretical curve agrees very well with experimental data below the glass transition 
temperature Tg (see the experimental data of Sosman |Q and Richet et al. |61[], and the 



theoretical curve of Horbach et al. in the figure). From the graph we see that an 
extrapolation of c^'' to lower temperatures extrapolates nicely to the experimental data and 
that an extrapolation of Cvib to lower temperatures can be joined smoothly to the curve from 
the harmonic approximation, thus showing that the two types of calculations are consistent 
with each other. 

The expressions derived in Sec. H were valid for all wave- vectors q and only at the end, 
i.e. in Eq. (|19|), we restricted ourselves to g = in order to obtain the equation relating the 
frequency dependent specific heat to the temperature fluctuations. After having investigated 
so far the temperature dependence of 0(2;), we now turn our attention to $tt(q', ^), the 
autocorrelation function of 6Tg{t), which measures the fluctuations in temperature at flnite 
wave- vectors. From the deflnition of ^TT{(l,t) it is clear that this function should scale like 
T^. That this is indeed the case is shown in Fig. |^, where we show ^TT{(l,t)/T^ for various 
wave-vectors and temperatures. Since for each wave-vector the curves for the different 
temperatures fall on top of each other we recognize immediately that the dependence 
is correct. Note that this weak temperature dependence for g > is in strong contrast to 
the one found for g = (see Fig. 0b). It reflects the fact that the thermal conductivity A is 
only a weak function of temperature, see Eq. (p6|). 

From the plot we also see that the typical time scale over which the correlation functions 
decay, increases with decreasing wave-vector, in qualitative agreement with Eq. (|2^) which 
predicts a dependence. To determine the g-dependence of this decay we deflne a decay 
time r(g) as the time it takes the correlation function to decay to 0.1 of their initial value. 
The wave- vector dependence of r(g) is shown in Fig. |^ where we plot g^T(g) versus g^. (Since 
within the accuracy of our data the temperature dependence of $tt(?5^) is independent of 
T we show only one set of data points.). From this flgure we recognize that for small wave- 
vectors the relaxation time scales indeed like g~^ + const, as expected from hydrodynamics 
(straight line), that however, this linear dependence breaks down for large wave- vectors. 
Furthermore we see that the slope of this straight line is positive, which means that the 
second term in Eq. (^) is larger than zero. From this equation it follows that the intercept 
of the straight line with the abscissa is given by nCp'^/X, where n is the particle density 
and A is the thermal conductivity. We read off an intercept 0.0180 ps/A^ and with the 
speciflc heat of 4kB per particle, see Fig. we obtain A = 2.4W/Km. This value is in 
good agreement to the one determined by a completely different method in the simulation 



by Jund and Jullien who found A ~ 1.3W/Km around lOOOK [Q. The experimental values 



for this quantity range between 2 and 3W/Km at high temperatures if T > lOOOi^' |^ 
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i.e. our value is in agreement also with the experimental data. (Note that in experiments 
it is found that the A is a strong function of temperature for temperatures below ^ lOOOK. 
For higher temperatures it seems, however, to level off and thus it can be extrapolated 
reasonably safely to temperatures in the melt. Since this temperature dependence of A is 
due to anharmonicities one can conclude that these become effectively independent of T in 
the temperature range of the supercooled melt.) 



V. SUMMARY AND DISCUSSION 

The goal of this paper is to show how the frequency dependent specific heat, Cv{z), is 
related to the dynamics of the particles on the microscopic level. For this we use the Mori- 
Zwanzig projection operator formalism to derive an exact expression for Cy{z) (Eq. (p!7D). 
This expression allows us to identify the physical mechanism which causes the frequency 
dependence of c^,(z), namely the relaxation of the potential energy during the structural re- 
laxation. Using an exact transformation formula by Lebowitz et ai, we obtain an equation 
which relates Cy{z) to the Laplace transform of K{t), the autocorrelation function of tem- 
perature fluctuations (Eq. (p^)) in the microcanonical ensemble, and which thus can be used 
to determine Cy{z) from a computer simulation. This relation has been derived previously 
by Nielsen on the basis of thermodynamic arguments. In contrast to that approach we 
are, however, also able to generalize the correlator K{t) to finite wave-vectors and to relate 
the time dependence of these quantities to the thermal conductivity of the system, Eq. (p6[). 

By using molecular dynamics computer simulations of a simple but quite realistic model 
for silica, we have determined the time and temperature dependence of K{t). We see that at 
low temperatures this correlator shows a two-step decay, similar to the one that is found in 
the time correlation functions for structural quantities, such as the intermediate scattering 
function. In contrast to these correlators the height of the plateau at intermediate times 
decreases with decreasing temperature, a trend that can be understood by realizing that at 
very low temperatures K{t) is directly related to the autocorrelation function of the velocity. 

From the time dependence of K{t) we have calculated the frequency dependent specific 
heat. In contrast to previous numerical investigations the accuracy of our data is high 
enough to analyze in detail the frequency dependence of the real and imaginary part of 
c{z). We find that at very high frequencies the value of c'{u) is the one of an ideal gas 
and that with decreasing frequency it shows a rapid increase which is due to the vibrational 
excitations of the system. At low temperatures we see that c'{lj) shows a second increase at 
frequencies which correspond to the time scales of the structural relaxation. This frequency 
dependence is also reflected in c"{uj) where the first and second increase are reflected by the 
microscopic and a— peak, respectively. Thus we find that the frequency dependence of c{z) 
is qualitatively very similar to the one found in the dynamical susceptibility for structural 
quantities, which shows how intimately connected these quantities are. Further evidence for 
this can be obtained from the observation that the location of the a— peak in c"{uj) shows the 
same temperature dependence as the structural relaxation, in agreement with experimental 



findings [10,12,13 



From the height of the two mentioned steps in c'{u), we are able to determine the 
vibrational and the configurational part of the specific heat. We find that the former is 
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significantly higher than the latter, which is in agreement with the experimental observation 
that in strong liquids the drop in the specific heat at the glass transition is relatively small. 

Finally we have calculated the time dependence of the autocorrelation functions of tem- 
perature fiuctuations at finite wave-vectors. In agreement with our theoretical prediction, 
these functions decay much faster than the one for q = 0, i.e. K{t), and depend only very 
weakly on temperature. From the g— dependence of the relaxation time of this correlator we 
calculate the thermal conductivity A and find good agreement of our value with the one in 
experiments and a computer simulation in which A was determined by a different method. 

We also point out that, since our simulations have been done at constant volume, it 
is clear that the frequency dependence of c{z) and the strong temperature dependence of 
Umax, the location of the a— peak in c"{u!), is noi the result of the frequency and temperature 
dependence of the macroscopic density. Some time ago Zwanzig proposed (essentially) the 



following mechanism for the T— dependence of cOmax ^ change in temperature will in 
general give rise to a change in density (since most real experiments are done at constant 
pressure and not density) . Due to the high value of the bulk viscosity, this volume relaxation 
will be slow and occur on the time scale of the a— relaxation, and hence Cp will be frequency 
dependent. Since in turn the frequency dependence of the viscosity is due to the slow 
relaxation of the structure on the microscopic scale, Zwanzig thus argued that the reason 
for the T— dependence of Umax is just an indirect effect of the slow microscopic dynamics. 
Since in a system with constant volume this mechanism is not present and our simulations 
demonstrate that Umax does show a strong T— dependence, we conclude that the reason for 
this dependence must be a different one. 

Finally we also mention that from the shape of the a— peak in c"{uj) it is also possible to 
calculate the time dependence of the enthalpy in a cooling and heating experiment. For this 
we made simulations in which the sample was first cooled with a finite cooling rate through 
the glass transition temperature and subsequently reheated above Tg. Using the equilibrium 
data for c{z) we were able to reproduce accurately the time and temperature dependence 
of the enthalpy , which shows that if one knows the equilibrium quantity c{z) , one is 

also able to predict the system in an out-of-equilibrium situation. 
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APPENDIX A: RELATION BETWEEN THE DENSITY OF STATES AND THE 
AUTOCORRELATION FUNCTION OF THE KINETIC ENERGY FOR A 

HARMONIC SYSTEM 

For a purely harmonic system with the Hamiltonian 

H = ^ + ^ (Al) 
2m 2 ^ ' 

momenta and space coordinates are Gaussian variables. This simplifies considerably the 
evaluation of the normalized autocorrelation function of the kinetic energy K{t). For a set 
of Gaussian variables with zero mean the four point correlation function (ABCD) can be 
expressed by the two point correlation functions: 
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{ABCD) = {AB) {CD) + {AC) {ED) + {AD) {BC) . (A2) 
Using this relation, the autocorrelation function {p^^(t)pl) can be written as 

{Plit)pl) = 2{p,m)' + {pl){pl) . (A3) 
Since {p^{t)py) = 5^u{p'j^ cos fit, the autocorrelation function $Tr(^) is given by 

1 



2 



Y.M)P,? (A5) 



— (cos(2fit) + 1) . (A6) 



All averages are in the canonical ensemble and we used that (p^) = ksT . From Eq. (^TD 
we know, by using the value of the specific heat of a harmonic oscillator in three dimensions 
d"^ = 3kB, that $jT^''(t) = $rr(t) — T^/S. By defining the velocity autocorrelation function 
J{t) = m(v(i(:)v) and noting, that for a harmonic oscillator {v{t)\v) = S/c^T/m cos(fit) we 
obtain the result 



Taking the Fourier transform of ( [A7| ) and using, that the density of states g\ 
2J{uj)/3kBT we arrive at the result 



to 



i^M = ^ . (A8) 
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FIG. 1. Time dependence of the normalized kinetic energy autocorrelation function, a) 
T = 3000K; each thin curve is the average over ten different samples; the average over these 
curves gives the bold curve, b) K{t)/K{{)) for all temperatures investigated, c) Bold curves: K{t) 
for T = 3000K and T = 300K; thin curves: velocity-autocorrelation function J{2t)/6kBT (see 
Appendix A) at the same temperatures. 
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FIG. 2. Frequency dependence of K{lo) and g{u)/2)/9> (dashed and solid curves, respectively). 
Bold lines: T = 3000K, thin lines: T = 300K. 
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FIG. 3. Frequency dependence of the specific heat for all temperatures investigated, a) Real 
part, b) Imaginary part. 



24 



0.5 




10"^ 10"^ lO""" 10"^ 10"^ 10"' 10° 

co/27r [THz] 

FIG. 4. Main figure: Frequency dependence of the a— peak for intermediate and low temper- 
atures. Inset: S ame curves scaled by the height of their maximum versus u> j Wmax ) where Wmax is 
the location of the maximum. Dashed line: Kohlrausch- Williams- Watts function. 
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FIG. 5. Main figure: Temperature dependence of the a— relaxation times as determined from 
the incoherent intermediate scattering function and the frequency dependent specific heat. The 
straight lines are fits with an Arrhenius law. Inset: ratio of these relaxation times. 
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FIG. 6. Temperature dependence of Cconf and Cvibi the configurational and vibrational part of 
the specific heat, respectively, and of their sum c^^ . The dashed line is the specific heat as calculated 
from the harmonic approximation |^3|. The curves with the small symbols are experimental data 
from Ref. 
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FIG. 7. Time dependence of the autocorrelation function of the generalized temperature fluc- 
tuations 6Tq{t). The different line styles correspond to different temperatures and the different 
thickness to different wave-vectors. 
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FIG. 8. Relaxation time for the autocorrelation function ^TT{Q,t)/T'^ multiplied by q'^ versus 
the square of the wave-vector. The straight line is a fit with the functional form given by Eq. (| 
This line is given by q'^r = O.OlSOpsA-^ + O.OlOlg^ps. 



27 



